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Abstract 

It has been observed that the earthquake events possess short-term memory, i.e. that events 
occurring in a particular location are dependent on the short history of that location. We con- 
duct an analysis to see whether real-time earthquake data also possess long-term memory and, if 
so, whether such autocorrelations depend on the size of earthquakes within close spatiotemporal 
proximity. We analyze the seismic waveform database recorded by 64 stations in Japan, includ- 
ing the 2011 "Great East Japan Earthquake", one of the five most powerful earthquakes ever 
recorded which resulted in a tsunami and devastating nuclear accidents. We explore the question 
of seismic memory through use of mean conditional intervals and detrended fiuctuation analysis 
(DFA). We find that the waveform sign series show long-range power-law anticorrelations while 
the interval series show long-range power-law correlations. We find size-dependence in earthquake 
auto-correlations — as earthquake size increases, both of these correlation behaviors strengthen. We 
also find that the DFA scaling exponent a has no dependence on earthquake hypocenter depth or 
epicentral distance. 

PACS numbers: PACS numbers:89.65.Gh, 89.20.-a, 02.50.Ey 
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I. INTRODUCTION 



Many complex physical systems exhibit complex dynamics in which subunits of the system 
interact at widely varying scales of time and space l|, |2| . These complex interactions often 
generate very noisy output signals which still exhibit scale-invariant structure. Such complex 
systems span areas studied in physiology js], finance j4|, and seismology js-li]. 

In seismology the study of seismic waves is both scientifically interesting and of practical 
concern, particularly in such applied areas as engineering. A better understanding of seis- 
mic waves is immediately applicable in the design of structures for earthquake-prone areas 
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17|. It also allows scientists to better understand the underlying mechanisms that drive 



earthquakes 



18|-|23[|. In seismology, temporal and spatial clustering are considered impor- 



tant properties of seismic occurrences and, together with the Omori law (dictating aftershock 
timing) and the Gutenberg- Richter law (specifying the distribution of earthquake size), com- 
prise the main starting requirements to be fulfilled in any reasonable seismic probabilistic 
model. Analyzing the timing of individual earthquakes, Ref. [sf introduces the scaling con- 
cept to statistical seismology. The recurrence times are defined as the time intervals between 
consecutive events, Tj = tj — In the case of stationary seismicity, the probability density 
P(r) of the occurrence times was found to follow a universal scaling law 



P(r) = RfiRr) 



(1) 



where / is a scaling function and R is the rate of seismic occurrence, defined as the mean 
number of events with M > Mc 0]. Reference 0, ^ has demonstrated how the structure of 
seismic occurrence in time and magnitude can be treated within the framework of critical 
phenomena. 

Recently, a few papers have analyzed the existence of correlations between magnitudes 
of subsequent earthquakes [3, B]. Analyzing earthquakes with r greater than 30 minutes, 
Ref. 71] reported possible magnitude correlations in the Southern California catalog. Magni- 
tude correlations have often been interpreted as a spurious effect due to so called short-term 

n 

aftershock incompleteness (STAI) [9|. This hypothesis assumes that some aftershocks, espe- 
cially small events, are not reported in the experimental catalogs, which is in agreement with 
the standard approach that assumes interdependence of earthquake magnitudes implying no 
memory in earthquakes. 
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However, recent work has also challenged this interpretation. Reference lOj] reports the 
existence of magnitude clustering in which earthquakes of a given magnitude are more likely 
to occur close in time and space to other events of similar magnitude. They find that a 
subsequent earthquake tends to have a magnitude similar to but smaller than the previous 
earthquake. Reference [81] also reports the existence of magnitude correlations and addi- 
tionally demonstrates the structure of these correlations and their relationship to At and 
Ar, where the latter represents the distance between subsequent epicenters. Reference I2I 
creates a model to explain these magnitude correlations. They note that the Omori law and 
"background tectonic cycles" are responsible for clustering in interoccurrence times. Addi- 
tionally, Refs. 13 1 and [l^ find that the distribution of recurrence times strongly depends 
on the previous recurrence time such that small and large recurrence times tend to cluster 
in time. This dependence on the past is reflected in both the conditional mean recurrence 
time and the conditional mean residual time until the next earthquake. 

Since it is our hypothesis that long-range autocorrelations exist in seismic waves, we first 
note that long-range power-law autocorrelations are quite common in a large number of 
natural phenomena ranging from weather 2j426|, and physiological systems [3|, l27H30|. to 



financial markets 



31 



-l37|. 



In addition to analyzing the raw waveform, it is also common to analyze related time 
series, such the time series generated by taking the sign or magnitude of the wave- 

reports an empirical approximate relation at small time scales 
for the scaling exponents calculated for sign, magnitude, and the original time series, 
ttsign = 1/2(0; magnitude + ^original), physiology. The study of magnitude and sign time 
series is important in physiology because the magnitude time series exhibits weaker auto- 
correlations and a scaling exponent closer to the exponent of an uncorrelated series found 
when a subject is unhealthy |3|. Diagnostic power in physiology has been confirmed for sign 
time series as well — the sign time series of heart failure subjects exhibit scaling behavior 
similar to that observed in the original time series, but significantly different that of healthy 
subjects j^. Understanding the correlation properties of these three time series allows us to 
also understand the underlying processes generating them. 

Our investigation and discussion is organized as follows. First, we study the autocorre- 
lations of interval series by using the mean conditional technique. Second, we employ de- 
trended fluctuation analysis (DFA) 38|-|40| and find long-range power-law autocorrelations 
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in the sign and interval time series. For the interval time series we find a positive regression 
between the DFA scaling exponent a and earthquake size (measured by the Richter mag- 
nitude scale M or seismic moment Mq), while for the sign time series we find an inverted 
regression between a and earthquake magnitude. Thus we report that the observed autocor- 
relation depends on earthquake size, both in the sign and interval time series. We also find 
that the scaling exponent a has no dependence on hypocenter depth or epicentral distance. 



II. DATA 



Seismic waves are unique in that they have non-stationarities of a much larger order than 
those of any other known natural signal. Large earthquakes are characterized by a maximum 
amplitude that is often > 100 times larger than the mean amplitude [see Fig. 1(a)]. This is a 
limitation that makes seismic waves difficult to analyze using traditional analysis. Although 
we might want to use detrended ffuctuation analysis (DFA) |38l-l41| . originally proposed 
to study the correlations in a time series in the presence of non-stationarities commonly 
observed in natural phenomena, the level of non-stationarity in earthquakes is so large that 



DFA is inappropriate regardless of the order of the polynomial fit applied [40|]. Thus, due 
to lack of methods for highly non-stationary signals, we do not analyze correlations in the 
series of magnitudes, but instead analyze the correlations in the sign series [Fig. 1(c)] and 
interval series [Fig. 1(d)]. For our data, we use the seismic waveform database from the 
National Research Institute for Earth Science and Disaster Prevention (NIED) F-net (Full 
Range Seismograph Network of Japan), which records continuous seismic waveform data Wt 
by using broadband sensors in 64 stations in Japan [see Fig. 1(a)]. In our study we select 
46 stations (ADM, AOG, ASI, HID, HJO, HRO, IGK, IMG, INN, lYG, IZH, KGM, KMU, 
KNM, KNP, KNY, KSK, KSN, KSR, KYK, MMA, NKG, NOK, NOP, NRW, NSK, OSW, 
SAG, SHR, SIB, TAS, TGA, TGW, TKO, TMC, TSA, TYM, TYS, UMJ, WTR, YAS, 
YNG, YSI , YTY, YZK, ZMM), based on locations and integrity of data series. Seismic 
signals are recorded in three directions: (1) U (up-down with up positive), N (north-south 
with north positive), and E (east- west with east positive) 20]. In this paper, we report 
results from the vertical dimension only (U data), since the results for the horizontal data 
(N and E) data are very similar. Sampling intervals have five recording frequencies: 80Hz, 
20Hz, IHz, O.lHz, and O.OlHz. We study earthquake coda wave data with IHz sampling 
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interval for the year 2003, together with selected earthquake coda wave data from 11 March 
2011. We note that, because of the interaction between earthquakes, not all earthquakes 
can be employed in our analysis (see Appendix A). The data from 11 March 2011 is selected 
because it contains the notable 2011 Tohoku earthquake ("Great East Japan Earthquake") 
which resulted in the tsunami that caused a number of nuclear accidents. We also add two 
large earthquakes (M = 7.3 and M = 7.6) to our study, which also occurred the same day 
as aftershocks. 

We employ the following procedure to create our time series: 

(i) For each selected earthquake (see Appendix A) we create a new time series, the nor- 
malized waveform denoted by Wt out of the raw seismic acceleration waveform data 



Wnorm = {Wt-w)/yw'^ -W"^. (2) 

(ii) From the time series Wnorm we define a new sub-series starting at time coordinate 
where maximum wt occurs and terminating at the end of the normalized waveform 
(see inset in Figl(a)). 

(iii) Let the time series tj denote the points in time when w'^ changes sign, with tj < tj+i. 
We define (see Fig 1(c)) the interval series by 

Ti = ti- U-i. (3) 

(iv) The sign series (see Fig. 1(d)) is defined by 

St = sgn{w[) (4) 

Note that our definition of interval is different than that recently defined in several papers, 
where the return intervals r have studied between consecutive fluctuations above a volatility 
threshold q in different complex systems. The probability density function (pdf) of return 
intervals Pqir) scales with the mean return interval as 

P,{r) = r-'fir/T) (5) 



where /() is a stretched exponential |32l434| . Since, on average, there is one volatility above 



the threshold q for every Tq volatilities, then it holds that [36 1 

/"OO 

1/r, ^ / Pi\R\)d\R\ = Pi\R\ > g) ~ g-". (6) 

Jq 



For the time intervals between events given by fluctuations R where R > q Ref. 36 1 
derived that r^, the average of Tq, obeys a scaling law, 

^ = (7) 

where by a denotes our estimate of the tail exponent probability density function, P(|i?|^+"). 
Similarly, if P{\R\) follows an exponential function P{\R\) oc exp(— then employing 
Eq. ([6]) we easily derive 

Tg OC exp (8) 
Eq. (IH]) can be used as a new method for estimation of the exponential parameter /3. 



III. MEMORY OF INTERVAL TIME SERIES 



Returning to waveform data, we begin analyzing the series by studying the conditional 
mean 

{r\ro)/T (9) 

which gives the mean value of r (see Eq. ([3])) immediately following a given term tq, normal- 
ized in units of r. The conditional mean gives evidence of whether seismic memory exists in 
the intervals in the form of correlations or anticorrelations. For example, should correlations 
exist, one would expect the mean interval to be shorter in the window immediately following 
a small interval. 

Indeed, Fig. 2 shows that the large intervals r tend to follow large initial tq and small r 
follow small tq indicating the existence of (positive) correlations in the interval time series. 
We also note that the autocorrelations tend to be stronger for the subset associated with 
larger earthquakes than for those associated with smaller earthquakes. 

To expand on this we also extend our investigation to longer range effects. We investigate 
the mean interval after a cluster of n consecutive intervals that are either entirely above the 
series mean or entirely below it. We denote clusters that are entirely above the series mean 
with a "+" and clusters below the series mean with a "— " . Fig. 3 shows the mean interval 
T that follows a To(n) defined as a cluster size of n. We find that for "+" clusters — shown 
by open symbols — the mean interval increases with the size of the cluster n. This is the 
opposite of what we find for "— " clusters — shown as closed symbols. The results indicate 
the existence of at least short-term memory in the interval time series. Furthermore, we 
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find that the mean interval increases with the seismic magnitude. However, this relationship 
breaks at the high end of the Richter magnitude scale M > 6.5. 



IV. DETRENDED FLUCTUATION ANALYSIS 



Many physical, physiological, biological, and social systems are characterized by complex 
interactions between a large number of individual components, which manifest in scale- 
invariant correlations [2I, Q, 43|. Since the resulting observable at each moment is the 
product of a magnitude and a sign, many recent investigations have focused on the study 



of correlations in magnitude and sign time series J3 
time series of changes Sti of heartbeat intervals [3 



intratrading times in the stock market 




44j-|46| . For example, the 



46| . physical activity levels 



44l ]. and river flux values 45| all exhibit power-law 



anticorrelations, while their magnitudes \6Ti\ are positively correlated. A common means 
of finding autocorrelations hidden within a noisy non-stationary time series is detrended 
fluctuation analysis (DFA) 38|-|40| . In the DFA method, the time series is partitioned into 



pieces of equal size n. For each piece, the local trend is subtracted and the resulting standard 
deviation over the entire series is obtained. In general, the standard deviation F{n) of the 
detrended fluctuations depends on n, with smaller n resulting in trends that more closely 
match the data. The dependence of F on can generally be represented as a power law 
such that 

F{n) oc n", (10) 

where a is the scaling exponent — sometimes referred to as the Hurst exponent — to be ob- 
tained empirically. DFA therefore can conceptually be understood as characterizing the 
motion of a random walker whose steps are given by the time series. F{n) gives the walker's 
deviation from the local trend as a function of the trend window. Because the root mean 
square displacement of a walker with no correlations between his steps scales like \/\n), we 
can expect a time series with no autocorrelations to yield an a of 0.5. Similarly, long-range 
power-law correlations in the signal (i.e. large terms follow large terms and small terms 
follow small terms) manifest as a > 0.5. Power-law anticorrelations within a signal will 
result in a < 0.5. Additionally, DFA can be related to the autocorrelation as follows: if the 
autocorrelation function C{L) can be approximated by a power law with exponent 7 such 
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that 

C{L) oc L-"', (11) 



then 7 is related to a by 



38| 

a ^1-7/2. (12) 



Another reason we employ the DFA method is that it is appropriate for sign time series 



281]. Other techniques for the detection of correlations in non-stationary time series are not 



appropriate for sign time series. Also, because the sign and interval time series have affine 
relations, the analysis of sign will be helpful in understanding the intervals. However, the 
DFA gives biased estimates for the power-law exponent in analysis of anticorrelated series 



391], and so in order to improve the accuracy of analysis, we integrate the time series before 



we employ the standard DFA procedure. 

For the 2011 Tohoku earthquake, also known as the "Great East Japan Earthquake", 
we present the fluctuation function F{n) of the coda wave, measured at KSN station, as 
typical examples of sign and intervals time series (Fig. 4). By using DFA, we find, for most 
coda waves after earthquakes, that the time series of the intervals are consistent with a 
power-law correlated behavior a = 0.69, while the sign time series of Eq. (jll) are consistent 
with a power-law anti-correlated behavior {a = 0.32). The results therefore indicate that for 
the interval series large increments are more likely to be followed by large increments and 
small increments by small increments. These results are in agreement with the results of 
the correlation analysis reported in Section 3. In contrast, anticorrelations in the sign time 
series indicate that positive increments are more likely to be followed by negative increments 
and vice versa. 

For the entire set of sign time series comprising our sample we calculate the average DFA 
scaling exponent a = 0.34 ± 0.09 indicating anticorrelations, and for the interval time series 
we calculate the average DFA scaling exponent a = 0.58 ± 0.08 indicating correlations. For 
the different stations measuring the 2011 Tohoku earthquake we find that for the sign time 
series, a = 0.29 ± 0.05 and for the interval time series, a = 0.66 ± 0.07. 
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V. RELATION BETWEEN EARTHQUAKE MOMENTS AND SCALING EXPO- 
NENTS OF SIGN AND INTERVAL SERIES 



Because large earthquake events release such extraordinary amounts of energy, it is rea- 
sonable to ask whether their occurrence influences local wave dynamics. To this end, we 
study interval time series of coda waves from earthquakes occurring in 2003, also including 
the particularly large events of 11 March 2011, when three events M > 7 occurred in the 
same day. Fig. 5(a) shows the DFA scaling exponent of the sign series versus seismic mo- 
ment, where seismic moment is a quantity used to measure the size of an earthquake. We 
flnd a decreasing functional dependence between the DFA exponent of the sign series and 
the seismic moment of the proximal earthquake with slope 7 = —0.028 ± 0.002, indicating 
that the DFA exponent decreases approximately with seismic moment. Note that because 
most of the exponents arc < 0.5, this indicates the presence of ever stronger anticorrclations 
in the time series as earthquake magnitude increases. Note, however, that the data break 
with this trend for very large earthquakes (Richter magnitude scale > 6.6 or seismic moment 
> 10^9). 

We also find similar results in the interval series, the difference being that the anticor- 
rclations become correlations. Fig. 5(b) shows that the DFA interval exponent and seismic 
moment exhibit a positive functional dependence with slope 7 = 0.025 ± 0.002 so that the 
DFA exponent increases with increasing seismic moment. Because most of the exponents 
for the interval series are > 0.5, this indicates that the series show stronger correlations for 
increasing seismic moment. Again, as with the sign series, we find a deviation from this 
trend for very large earthquakes. 

Having observed the influence of seismic moment on autocorrelations, we now investigate 
whether other readily observable factors such hypoccntcr depth and cpiccntral distance (the 
distance from the event to the recording station) also contribute. Specifically, we would like 
to explore whether there is evidence that such long-term memory is affected by the spreading 
process as seismic waves disseminate outward from their epicenter to a recording station or 
whether the memory observed is strictly due to the seismic activity. Fig. 6 shows that the 
DFA exponent for both interval and sign series are independent of both hypocenter depth 
and epicentral distance. From these results we speculate that the DFA exponent is mainly a 
result of the characteristics of the hypocenter rather than the process by which the seismic 
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waves are spread. 

For moderately large earthquakes (Mq = 10^^ ~ 10^^), we approximate the relation 
between the DFA scaling exponent and seismic moment through the empirical formula 

a^a logio{Mo) + c (13) 

where a = —0.028, c = 0.797 for the sign time series and where a = 0.025, c = 0.174 for 
the interval time series. Since 

M=(/o(?(Mo)-9.1)/1.5, (14) 

we can also write 

a ^ a{1.5M + 9.1) + c = a' M + c', (15) 

where a' = —0.042, c' = 0.542 for the sign series, and a' = 0.037, c' = 0.398 for the 
interval series. 



We note that similar size dependence in Hurst exponent was found in Ref. 47|] where 



Hurst exponents of financial time series increase logarithmically with company size. 
VI. SUMMARY 

We analyze seismic coda waves during earthquakes, finding long-range power-law auto- 
correlations in both the interval and sign time series. The sign series generally display power- 
law anticorrelated behavior, with anticorrelations becoming stronger with larger earthquake 
events, while the interval series generally display power-law correlated behavior, with corre- 
lations also becoming stronger with larger earthquake events. We also show that while the 
DFA autocorrelation exponent is influenced by the size of the earthquake seismic moment, 
it is unaffected by earthquake depth or epicentral distance. Our findings are in contrast 
with a standard approach which assumes independence in earthquake signals and thus have 



strong implications on the ongoing debate about earthquake predictability 481]. 
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VIII. APPENDIX: THE SELECTION OF EARTHQUAKES 

In some regions it is common for multiple earthquakes to occur in short succession. In 
many cases, because the interoccurrence times are so short, the coda waves can be derived 
from more than one earthquake. This is especially true for large earthquakes with many 
aftershocks 49(|. In order to make sure that the coda waves we study are the effects of 
only one earthquake, we need a way of determining which earthquakes are independent. We 
use the following two functions to determine the sphere of influence and duration of each 



earthquake by using the Richter magnitude scale M [49|. We select only those earthquakes 



that have no larger earthquake in their spatiotemporal sphere of influence, 

t ^ 10(A^-4-71)/l-67 (16) 

and 

R^2x 10(^^+i)/2-7, (17) 

where t is the duration and R is the sphere radius of influence. The two functions are 
empirical formulas based on an analysis of earthquakes in Japan \^\. The 10^^^^ '^/2.7 is an 
empirical formula that indicates the maximum radius that a human can feel an earthquake, 
especially for the earthquakes in Japan. 
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FIG. 1: (a) Location map for the 46 broadband stations of Full Range Seismograph Network of 
Japan (F-net) (red snow marks). Inset: An example of a record of a seismic wave (Up-Down 
component), (b) A part of the coda wave series indicated in inset of (a), as an example, (c) 
An example sign time series where the positive 4ign (+1) represents a positive waveform, and the 
negative sign (-1) represents a negative waveform in coda wave series of seismic wave, (d) Interval 
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FIG. 2: Scaled mean conditional interval {t\tq) /t vs tq/t . Five groups, one with no proximal 
earthquake and earthquakes with Richter magnitude scale M < 4.5, M = 4.5 ~ 5.5, M = 6.5 ~ 6.5, 
M > 6.5. An increasing trend implies a short-range correlation in the interval series. 
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FIG. 3: Long-range memory in interval clusters, tq signifies a cluster of intervals, consisting of n 
consecutive values that all are above (denote as "+") or below (denote as "— ") the median of the 
entire interval records. Plots display the scaled mean interval conditioned on a cluster, {t\tq) /t 
vs the size n of the cluster for five group intervals. The upper part (overplotted) of curves is 
for "+" clusters while the lower part is for "— " clusters. The plots show that "+" clusters are 
likely to be followed by large intervals and "— " clusters by small intervals, consistent with long- 
term correlations in interval records. Similar to Fig. 2, the long-term correlation increases with 
earthquake size, with exceptions for very large earthquakes. 
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FIG. 4: DFA fluctuation function F{n) of 2011 Tohoku earthquake as a function of time scale n 
{F{n) oc n") for (a) sign time series (a + 1 = 1.32, (a < 0.5), indicates anticorrelations) and (b) 
interval time series (a + 1 = 1.69, (a > 0.5), indicates correlations). 
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FIG. 5: Scaling exponent a vs seismic moment (Richter magnitude scale) for (a) sign time series 
(correlation coefficient Cor = —0.3604), and (b) interval time series (correlation coefficient Cor = 
0.3602). The values of 7 show negative slope in the regression a vs seismic moment of the sign 
series, and positive slope in the regression of the interval series. Triangular symbols show the mean 
of exponent within each bin ( bins: < le + 15, le + 15 ~ le + 16, le + 16 ~ le + 17, le + 17 ~ 
le+18, le+18 ~ le+19, le+19 ~ le+20, > le+21), the error bar shows the ± standard deviation. 
The plots show a linear relationship between logarithmic earthquake moment and scaling exponent 
a in the sign and interval series, with exceptions for very large earthquakes. 
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FIG. 6: Scaling exponent a vs hypocenter depth for events where Richter magnitude scale M < 5 
for (a) sign time series (b) interval time series. Inset: scaling exponent a vs hypocenter depth 
requiring that Richter magnitude scale M > 5. (c) and (d) show Scaling exponent a vs epicentral 
distance for events where Richter magnitude scale M < 5. Inset: scaling exponent a vs epicentral 
distance requiring that Richter magnitude scale M > 5. (c) sign time series, (d) interval series. 
All absolute values of correlation coefficient are smaller than 0.1, showing that a is uncorrelated 
with both hypocenter depth and epicentral distance. 
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